// Copyright (C) 2011  Andrew H. Chan
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

#include <iostream>
#include <iomanip>
#include <vector>
#include <cstdlib>
#include <set>
#include <algorithm>
#include <assert.h>

#include "ntSeq.hpp"
#include "differenceCoverSample.hpp"
#include "findMedian.hpp"

#include "blockBWT.hpp"

typedef unsigned int uint;
typedef unsigned long long uint64;

using namespace std;

BWT::BWT(ntSeq const &T, uint64 logv, uint64 numBlocks) : n(T.size()) {
    assert(T.size()>0);
    bwt.resize(n+1);

    cout << "numBlocks: " << numBlocks << endl;
    if(numBlocks != min(numBlocks,n)) {
        numBlocks = min(numBlocks,n);
        cout << "Reducing number of blocks to " << numBlocks << " because input string is of length " << n << endl;
    }

    for(uint64 ii=0; ii<4; ++ii) C[ii]=0;
    for(uint64 ii=0; ii<4; ++ii) {
        occ[ii] = vector<uint64>(n/INDEX_PERIOD+1);
        occ[ii][0] = 0;
    }
    pos = vector<uint64>(n/INDEX_PERIOD+1);

    //cout << "creating DCS" << endl;
    DCS dcs(T, logv);
    assert(dcs.n==n);

    // Generate splitters
    srand(0);
    std::vector<uint64> splitters = generateSplitters(T, dcs, numBlocks);
    cout << "numsplitters: " << splitters.size() << endl;

    uint64 charCounters[4]={0,0,0,0};

    // assume T is not empty string
    // handle suffix n (empty suffix)
    uchar newBase = T.get(n-1);
    bwt.set(0, newBase); 
    ++C[newBase];
    ++charCounters[newBase];
    pos[0]=n-1;

    uint64 bwtIndex=1;
    if(numBlocks>1) {
        // handle initial case
        uint64 firstSplitter=splitters[0];
        vector<uint64> smaller = smallerSuffixes(T, dcs, firstSplitter);
        assert(!verifySmallerSuffixes(T, firstSplitter, smaller));
        multiQsort(T, dcs, smaller, 0, 0, smaller.size());
        for(uint64 jj=0; jj<smaller.size(); ++jj) {
            if(bwtIndex%INDEX_PERIOD==0) {
                for(uint64 kk=0; kk<4; ++kk) {
                    occ[kk][bwtIndex/INDEX_PERIOD]=charCounters[kk];
                }
                pos[bwtIndex/INDEX_PERIOD] = smaller[jj]-1;
            }
            if(smaller[jj] == 0) {
                term = bwtIndex;
                uchar newBase = 0;
                bwt.set(bwtIndex,0);
            } else {
                uchar newBase = T.get(smaller[jj]-1);
                bwt.set(bwtIndex,newBase);
                ++C[newBase];
                ++charCounters[newBase];
            }
            ++bwtIndex;
        }

        for(uint ii=1; ii<splitters.size(); ++ii) {
            uint64 splitter1=splitters[ii-1];
            uint64 splitter2=splitters[ii];
            vector<uint64> betweenSuf = betweenSuffixes(T, dcs, splitter1, splitter2);
            assert(!verifyBetweenSuffixes(T, splitter1, splitter2, betweenSuf));

            multiQsort(T, dcs, betweenSuf, 0, 0, betweenSuf.size());
            assert(!verifyQsort(T, betweenSuf));
            for(uint64 jj=0; jj<betweenSuf.size(); ++jj) {
                if(bwtIndex%INDEX_PERIOD==0) {
                    for(uint64 kk=0; kk<4; ++kk) {
                        occ[kk][bwtIndex/INDEX_PERIOD]=charCounters[kk];
                    }
                    pos[bwtIndex/INDEX_PERIOD] = betweenSuf[jj]-1;
                }
                if(betweenSuf[jj] == 0) {
                    term = bwtIndex;
                    bwt.set(bwtIndex,0);
                } else {
                    uchar newBase = T.get(betweenSuf[jj]-1);
                    bwt.set(bwtIndex,newBase);
                    ++C[newBase];
                    ++charCounters[newBase];
                }
                ++bwtIndex;
            }
        }

        // handle final case
        uint64 lastSplitter=splitters[splitters.size()-1];
        vector<uint64> larger = largerEqSuffixes(T, dcs, lastSplitter);
        assert(!verifyLargerEqSuffixes(T, lastSplitter, larger));

        multiQsort(T, dcs, larger, 0, 0, larger.size());
        for(uint64 jj=0; jj<larger.size(); ++jj) {
            if(bwtIndex%INDEX_PERIOD==0) {
                for(uint64 kk=0; kk<4; ++kk) {
                    occ[kk][bwtIndex/INDEX_PERIOD]=charCounters[kk];
                }
                pos[bwtIndex/INDEX_PERIOD] = larger[jj]-1;
            }
            if(larger[jj] == 0) {
                term = bwtIndex;
                bwt.set(bwtIndex,0);
            } else {
                uchar newBase = T.get(larger[jj]-1);
                bwt.set(bwtIndex,newBase);
                ++C[newBase];
                ++charCounters[newBase];
            }
            ++bwtIndex;
        }
    } else { // special case for 1 block
        vector<uint64> allSuffixes(T.size());
        for(uint64 ii=0; ii<T.size(); ++ii) allSuffixes[ii]=ii;
        multiQsort(T, dcs, allSuffixes, 0, 0, allSuffixes.size());
        for(uint64 jj=0; jj<allSuffixes.size(); ++jj) {
            if(bwtIndex%INDEX_PERIOD==0) {
                for(uint64 kk=0; kk<4; ++kk) {
                    occ[kk][bwtIndex/INDEX_PERIOD]=charCounters[kk];
                }
                pos[bwtIndex/INDEX_PERIOD] = allSuffixes[jj]-1;
            }
            if(allSuffixes[jj] == 0) {
                term = bwtIndex;
                uchar newBase = 0;
                bwt.set(bwtIndex,0);
            } else {
                uchar newBase = T.get(allSuffixes[jj]-1);
                bwt.set(bwtIndex,newBase);
                ++C[newBase];
                ++charCounters[newBase];
            }
            ++bwtIndex;
        }
    }

    assert(bwtIndex==n+1);

    uint64 charSum=1;
    for(uint64 ii=0; ii<4; ++ii) {
        uint64 tmp=charSum;
        charSum+=C[ii];
        C[ii]=tmp;
    }

    cout << "bwt: ";
    for(uint64 ii=0; ii<bwt.size(); ++ii) {
        if(ii==term) cout << "$";
        else cout << ntSeq::invBM[bwt.get(ii)];
    }
    cout << endl;

    vector<uchar> unpermutedSeq = unpermute();
    assert(unpermutedSeq.size() == n);
    string unpermutedStr;
    for(uint64 ii=0; ii<unpermutedSeq.size(); ++ii)
        unpermutedStr+=ntSeq::invBM[unpermutedSeq[ii]];
    cout << "unpermuted: " << unpermutedStr << endl;
    assert(unpermutedStr == T.getSuffix(0));
}

vector<uint64> BWT::generateSplitters(ntSeq const &T, DCS const &dcs, uint64 numBlocks) const {
    assert(numBlocks>=1);
    set<uint64> tmpSplitters; // use set to avoid duplicate splitters
    while(tmpSplitters.size()<numBlocks-1) tmpSplitters.insert(rand()%n); // Using rand()
    vector<uint64> splitters(numBlocks-1);
    uint64 ii=0;
    for(set<uint64>::const_iterator citer=tmpSplitters.begin(); citer != tmpSplitters.end(); ++citer)
        splitters[ii++]=*citer;
    tmpSplitters.clear();
    // sort the splitters
    multiQsort(T, dcs, splitters, 0, 0, splitters.size());
    assert(!verifyQsort(T, splitters));
    return splitters;
}

vector<uint64> BWT::precomputePP(ntSeq const &T, uint64 v, uint64 pp) const {
    assert(pp<T.size());
    uint64 m=T.size()-pp;
    uint64 n=m-1;
    uint64 lenPP = min(n, v);

    vector<uint64> PP(lenPP);
    if(lenPP==0) {
        return PP;
    }
    uint64 kk=0, ll=0;
    // handle first case (i.e. lcp(P[0,m],P[1,m]))
    while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(pp+kk+1)) {
        ++kk;
        ++ll;
    }
    PP[0] = ll; // lcp(P[0,m],P[1,m])
    uint64 ii=1, jj=0;
    while(ii<lenPP) { // handle the rest
        if(ii>=kk) {
            kk=ii;
            ll=0;
        } else {
            ll=PP[ii-jj-1];
        }
        if(ii+ll==kk) {
            while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(pp+kk+1)) {
                ++kk;
                ++ll;
            }
            jj=ii;
        } else if(ii+ll>kk) {
            ll=kk-ii;
            jj=ii;
        }
        PP[ii] = ll; // P[ii] = lcp(P,P[ii+1,m])
        ++ii;
    }
    return PP;
}

vector<uint64> BWT::smallerSuffixes(ntSeq const &T, DCS const &dcs, uint64 pp) const {
    assert(T.size()>0);
    assert(pp<T.size());

    uint64 v = dcs.v;

    // precompute PP
    vector<uint64> PP = precomputePP(T, v, pp);
    uint64 n = T.size();
    uint64 m = n-pp;

    vector<uint64> smaller; // holds suffixes that are less than P
    uint64 kk=0, ll=0;
    while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) { // handle case lcp(T[0,n],P[0,m])
        ++kk;
        ++ll;
    }
    if(ll<v) { // compare mismatching characters
        if(ll!=m && (ll==n || T.get(ll)<T.get(pp+ll))) { // T[i..] < P
            smaller.push_back(0); // without reserving space in the vector, this could be slow
            assert(T.getSuffix(0) < T.getSuffix(pp));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(pp));
        }
    } else { // use difference cover sample
        if(dcs.lt(0,pp)) {
            smaller.push_back(0); // without reserving space in the vector, this could be slow
            assert(T.getSuffix(0) < T.getSuffix(pp));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(pp));
        }
    }

    uint64 ii=1, jj=0;
    while(ii<n) {
        if(ii>=kk) {
            kk=ii;
            ll=0;
        } else {
            assert(ii-jj-1<PP.size());
            ll=PP[ii-jj-1]; // PP[ii-jj] = lcp(P,P[ii-jj+1,m])
        }
        if(ii+ll==kk) {
            while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) {
                ++kk;
                ++ll;
            }
            jj=ii;
        } else if(ii+ll>kk) {
            if(kk-jj==v) {
                ll=kk-ii;
                while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) {
                    ++kk;
                    ++ll;
                }
                jj=ii;
            } else {
                ll=kk-ii;
                jj=ii;
            }
        }
        assert(ll<=v);
        if(ll<v) { // compare mismatching characters
            if(ll!=m && (ii+ll==n || T.get(ii+ll)<T.get(pp+ll))) { // T[i..] < P
                smaller.push_back(ii); // without reserving space in the vector, this could be slow
                assert(T.getSuffix(ii) < T.getSuffix(pp));
            } else {
                assert(T.getSuffix(ii) >= T.getSuffix(pp));
            }
        } else { // use difference cover sample
            if(dcs.lt(ii,pp)) {
                smaller.push_back(ii);
                assert(T.getSuffix(ii) < T.getSuffix(pp));
            } else {
                assert(T.getSuffix(ii) >= T.getSuffix(pp));
            }
        }
        ++ii;
    }

    return smaller;
}

int BWT::verifySmallerSuffixes(ntSeq const& T, uint64 pp, vector<uint64> const& smaller) const {
    cout << "Verifying smaller suffixes..." << endl;
    vector<uint64> suffixes(T.size(), false);
    // check that every suffix in smaller is less than pp
    for(uint64 ii=0; ii<smaller.size(); ++ii) {
        suffixes[smaller[ii]] = true;
        if(T.getSuffix(smaller[ii])>=T.getSuffix(pp))
            return 1;
    }
    // check that every suffix not in smaller is greater than or equal to pp
    for(uint64 ii=0; ii<T.size(); ++ii)
        if(!suffixes[ii])
            if(T.getSuffix(ii)<T.getSuffix(pp))
                return 1;
    return 0;
}

vector<uint64> BWT::largerEqSuffixes(ntSeq const &T, DCS const &dcs, uint64 pp) const {
    assert(T.size()>0);
    assert(pp<T.size());

    uint64 v = dcs.v;

    // precompute PP
    vector<uint64> PP = precomputePP(T, v, pp);
    uint64 n = T.size();
    uint64 m = n-pp;

    vector<uint64> larger; // holds suffixes that are less than P
    uint64 kk=0, ll=0;
    while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) { // handle case lcp(T[0,n],P[0,m])
        ++kk;
        ++ll;
    }
    if(ll<v) { // compare mismatching characters
        if(ll!=m && (ll==n || T.get(ll)<T.get(pp+ll))) { // T[i..] < P
            assert(T.getSuffix(0) < T.getSuffix(pp));
        } else {
            larger.push_back(0); // without reserving space in the vector, this could be slow
            assert(T.getSuffix(0) >= T.getSuffix(pp));
        }
    } else { // use difference cover sample
        if(dcs.lt(0,pp)) {
            assert(T.getSuffix(0) < T.getSuffix(pp));
        } else {
            larger.push_back(0); // without reserving space in the vector, this could be slow
            assert(T.getSuffix(0) >= T.getSuffix(pp));
        }
    }

    uint64 ii=1, jj=0;
    while(ii<n) {
        if(ii>=kk) {
            kk=ii;
            ll=0;
        } else {
            assert(ii-jj-1<PP.size());
            ll=PP[ii-jj-1]; // PP[ii-jj] = lcp(P,P[ii-jj+1,m])
        }
        if(ii+ll==kk) {
            while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) {
                ++kk;
                ++ll;
            }
            jj=ii;
        } else if(ii+ll>kk) {
            if(kk-jj==v) {
                ll=kk-ii;
                while(ll<m && ll<v && kk<n && T.get(pp+ll)==T.get(kk)) {
                    ++kk;
                    ++ll;
                }
                jj=ii;
            } else {
                ll=kk-ii;
                jj=ii;
            }
        }
        assert(ll<=v);
        if(ll<v) { // compare mismatching characters
            if(ll!=m && (ii+ll==n || T.get(ii+ll)<T.get(pp+ll))) { // T[i..] < P
                assert(T.getSuffix(ii) < T.getSuffix(pp));
            } else {
                larger.push_back(ii); // without reserving space in the vector, this could be slow
                assert(T.getSuffix(ii) >= T.getSuffix(pp));
            }
        } else { // use difference cover sample
            if(dcs.lt(ii,pp)) {
                assert(T.getSuffix(ii) < T.getSuffix(pp));
            } else {
                larger.push_back(ii);
                assert(T.getSuffix(ii) >= T.getSuffix(pp));
            }
        }
        ++ii;
    }

    return larger;
}

int BWT::verifyLargerEqSuffixes(ntSeq const& T, uint64 pp, vector<uint64> const& larger) const {
    cout << "Verifying larger suffixes..." << endl;
    vector<uint64> suffixes(T.size(), false);
    // check that every suffix in larger is larger than or equal to pp
    for(uint64 ii=0; ii<larger.size(); ++ii) {
        suffixes[larger[ii]] = true;
        if(T.getSuffix(larger[ii])<T.getSuffix(pp))
            return 1;
    }
    // check that every suffix not in larger is less than to pp
    for(uint64 ii=0; ii<T.size(); ++ii)
        if(!suffixes[ii])
            if(T.getSuffix(ii)>=T.getSuffix(pp))
                return 1;
    return 0;
}

vector<uint64> BWT::betweenSuffixes(ntSeq const &T, DCS const &dcs, uint64 splitter1, uint64 splitter2) const {
    vector<uint64> debugSmaller1;

    assert(T.size()>0);
    assert(splitter1<T.size());
    assert(splitter2<T.size());

    uint64 v = dcs.v;

    // precompute PP
    vector<uint64> PP1 = precomputePP(T, v, splitter1);
    vector<uint64> PP2 = precomputePP(T, v, splitter2);
    
    uint64 n = T.size();
    uint64 m1 = n-splitter1;
    uint64 m2 = n-splitter2;

    vector<uint64> betweenSuf; // holds suffixes that are greater than or equal to splitter 1 and less than splitter 2 
    uint64 kk1=0, ll1=0;
    while(ll1<m1 && ll1<v&& kk1<n && T.get(splitter1+ll1)==T.get(kk1)) { // handle case lcp(T[0,n],P1[0,m]) for splitter 1
        ++kk1;
        ++ll1;
    }
    uint64 kk2=0, ll2=0;
    while(ll2<m2 && ll2<v&& kk2<n && T.get(splitter2+ll2)==T.get(kk2)) { // handle case lcp(T[0,n],P2[0,m1]) for splitter 2
        ++kk2;
        ++ll2;
    }
    bool smaller1=false, smaller2=false;
    if(ll1<v) // compare mismatching characters
        if(ll1!=m1 && (ll1==n || T.get(ll1)<T.get(splitter1+ll1))) {
            smaller1=true; // T[i..] < P1
            assert(T.getSuffix(0) < T.getSuffix(splitter1));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(splitter1));
        }
    else // use difference cover sample
        if(dcs.lt(0,splitter1)) { 
            smaller1=true;
            assert(T.getSuffix(0) < T.getSuffix(splitter1));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(splitter1));
        }

    if(ll2<v) // compare mismatching characters
        if(ll2!=m2 && (ll2==n || T.get(ll2)<T.get(splitter2+ll2))) {
            smaller2=true; // T[i..] < P2
            assert(T.getSuffix(0) < T.getSuffix(splitter2));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(splitter2));
        }
    else // use difference cover sample
        if(dcs.lt(0,splitter2)) {
            smaller2=true;
            assert(T.getSuffix(0) < T.getSuffix(splitter2));
        } else {
            assert(T.getSuffix(0) >= T.getSuffix(splitter2));
        }

    if(!smaller1 && smaller2) betweenSuf.push_back(0); // Without reserving capacity in the vector beforehand, this could be slow 

    uint64 ii=1, jj1=0, jj2=0;
    while(ii<n) {
        // compute lcp for splitter 1
        if(ii>=kk1) {
            kk1=ii;
            ll1=0;
        } else {
            assert(ii-jj1-1<PP1.size());
            ll1=PP1[ii-jj1-1]; // PP1[ii-jj1] = lcp(P1,P1[ii-jj1+1,m1])
        }
        if(ii+ll1==kk1) {
            while(ll1<m1 && ll1<v && kk1<n && T.get(splitter1+ll1)==T.get(kk1)) {
                ++kk1;
                ++ll1;
            }
            jj1=ii;
        } else if(ii+ll1>kk1) {
            if(kk1-jj1==v) {
                ll1=kk1-ii;
                while(ll1<m1 && ll1<v && kk1<n && T.get(splitter1+ll1)==T.get(kk1)) {
                    ++kk1;
                    ++ll1;
                }
                jj1=ii;
            } else {
                ll1=kk1-ii;
                jj1=ii;
            }
        }

        // compute lcp for splitter 2
        if(ii>=kk2) {
            kk2=ii;
            ll2=0;
        } else {
            assert(ii-jj2-1<PP2.size());
            ll2=PP2[ii-jj2-1]; // PP2[ii-jj1] = lcp(P2,P2[ii-jj2+1,m2])
        }
        if(ii+ll2==kk2) {
            while(ll2<m2 && ll2<v && kk2<n && T.get(splitter2+ll2)==T.get(kk2)) {
                ++kk2;
                ++ll2;
            }
            jj2=ii;
        } else if(ii+ll2>kk2) {
            if(kk2-jj2==v) {
                ll2=kk2-ii;
                while(ll2<m2 && ll2<v && kk2<n && T.get(splitter2+ll2)==T.get(kk2)) {
                    ++kk2;
                    ++ll2;
                }
                jj2=ii;
            } else {
                ll2=kk2-ii;
                jj2=ii;
            }
        }

        assert(ll1<=v);
        assert(ll2<=v);
        bool smaller1=false, smaller2=false;
        if(ll1<v) { // compare mismatching characters
            if(ll1!=m1 && (ii+ll1==n || T.get(ii+ll1)<T.get(splitter1+ll1))) {
                smaller1=true; // T[i..] < P
                assert(T.getSuffix(ii)<T.getSuffix(splitter1));
            } else {
                assert(T.getSuffix(ii)>=T.getSuffix(splitter1));
            }
        } else {
            // use difference cover sample
            if(dcs.lt(ii,splitter1)) {
                smaller1=true;
                assert(T.getSuffix(ii)<T.getSuffix(splitter1));
            } else {
                assert(T.getSuffix(ii)>=T.getSuffix(splitter1));
            }
        }

        if(ll2<v) { // compare mismatching characters
            if(ll2!=m2 && (ii+ll2==n || T.get(ii+ll2)<T.get(splitter2+ll2))) {
                smaller2=true; // T[i..] < P
                assert(T.getSuffix(ii)<T.getSuffix(splitter2));
            } else {
                assert(T.getSuffix(ii)>=T.getSuffix(splitter2));
            }
        } else { // use difference cover sample
            if(dcs.lt(ii,splitter2)) {
                smaller2=true;
                assert(T.getSuffix(ii)<T.getSuffix(splitter2));
            } else {
                assert(T.getSuffix(ii)>=T.getSuffix(splitter2));
            }
        }
        if(!smaller1 && smaller2) betweenSuf.push_back(ii); // Without reserving capacity in the vector beforehand, this could be slow
        ++ii;
    }
    return betweenSuf;
}

int BWT::verifyBetweenSuffixes(ntSeq const& T, uint64 splitter1, uint64 splitter2, vector<uint64> const& betweenSuf) const {
    cout << "Verifying between suffixes..." << endl;
    vector<bool> suffixes(T.size(), false);
    for(uint64 ii=0; ii<betweenSuf.size(); ++ii) {
        suffixes[betweenSuf[ii]] = true;
        if(T.getSuffix(betweenSuf[ii])<T.getSuffix(splitter1) || T.getSuffix(betweenSuf[ii])>=T.getSuffix(splitter2)) {
            return 1;
        }
    }
    for(uint64 ii=0; ii<T.size(); ++ii) {
        if(!suffixes[ii]) {
            if(T.getSuffix(ii)>=T.getSuffix(splitter1) && T.getSuffix(ii)<T.getSuffix(splitter2)) {
                return 1;
            }
        }
    }
    return 0;
}

// v-order the suffixes in sufffixes
// modifies l
void BWT::multiQsort(ntSeq const &T, DCS const &dcs, vector<uint64> &suffixes, uint64 charIndex, uint64 index, uint64 length) const {
    if(length==0) return;

    uint64 charIncr = 21*charIndex;
    // selection sort serves as the base case
    if(length < 7) { // selection sort if partition has fewer than 7 elements 
        uint64 pLow=index, pLast=index+length-1;
        while(pLow < pLast) {
            uint64 minKey = dcs.packedKey(T,suffixes[pLow]+charIncr);
            uint64 pEqual = pLow+1;
            for(uint64 ii=pLow+1; ii<=pLast; ++ii) { // search for minimum element
                uint64 tmpKey = dcs.packedKey(T,suffixes[ii]+charIncr);
                if(tmpKey < minKey) {
                    minKey = tmpKey;
                    swap(suffixes[ii], suffixes[pLow]); // swap to beginning of selection group 
                    pEqual=pLow+1; // new min found, so reset pEqual to point to one past end of min elements
                } else if(tmpKey == minKey) {
                    swap(suffixes[ii], suffixes[pEqual]); // swap to min elements group at beginning of selection group
                    ++pEqual;
                }
            }
            uint64 numEqual = pEqual-pLow;
            
            if(numEqual>1) {
               if(charIndex<(dcs.v-1)/21) multiQsort(T, dcs, suffixes, charIndex+1, pLow, numEqual); // recursively sort equal group
               else qsortDCS(T, dcs, suffixes, pLow, numEqual);
            }
            pLow = pEqual; // prepare to sort next selection group
        }
        return;
    } else { // quicksort if partition has 7 or greater elements
        // choose pivot
        uint64 pivot;
        if(length <= 50) { // if partition has 50 or fewer elements, select the median of first, middle and last elements
            uint64 a = dcs.packedKey(T,suffixes[index]+charIncr), b = dcs.packedKey(T,suffixes[index+(length>>1)]+charIncr), c = dcs.packedKey(T,suffixes[index+length-1]+charIncr);
            pivot = a > b ? (b > c ? b : (a > c ? c : a)) : (a > c ? a : (b > c ? c : b));
        } else { // if partition has more than 50 elements, select median of 5 evenly spaced elements
            uint64 incr = length>>2;
            pivot = findMedian5(dcs.packedKey(T,suffixes[index]), dcs.packedKey(T,suffixes[index+incr]+charIncr), dcs.packedKey(T,suffixes[index+incr+incr]+charIncr), dcs.packedKey(T,suffixes[index+incr+incr+incr]+charIncr), dcs.packedKey(T,suffixes[index+length-1]+charIncr));
        }
        uint64 pCurLow = index, pLow = index, pCurHigh = index+length-1, pHigh = index+length-1;
        uint64 curKey;
        while(true) {
            while(pCurLow <= pCurHigh && (curKey=dcs.packedKey(T,suffixes[pCurLow]+charIncr)) <= pivot) {
                if(curKey == pivot) swap(suffixes[pCurLow],suffixes[pLow++]); // swap equal elements to beginning
                ++pCurLow;
            }
            while(pCurHigh > pCurLow && (curKey=dcs.packedKey(T,suffixes[pCurHigh]+charIncr)) >= pivot) {
                if(curKey == pivot) swap(suffixes[pCurHigh],suffixes[pHigh--]); // swap equal elements to end
                --pCurHigh;
            }
            if(pCurLow>=pCurHigh) break;
            swap(suffixes[pCurLow++],suffixes[pCurHigh--]); // swap low with high
        }
        pCurHigh = pCurLow-1;
        // pCurLow points one beyond the smaller elements
        // pCurHigh points one to the left of the larger elements

        // swap equal elements from left toward center
        uint64 numEqLeft = pLow-index; // number of elements equal to pivot
        uint64 numLow = pCurLow-pLow; // number of elements less than pivot
        if(numEqLeft > numLow)
            for(uint64 ii=index, jj=pLow, numSwap=numLow; numSwap>0; --numSwap, ++ii, ++jj) swap(suffixes[ii], suffixes[jj]);
        else
            for(uint64 ii=index, jj=pCurLow-numEqLeft, numSwap=numEqLeft; numSwap>0; --numSwap, ++ii, ++jj) swap(suffixes[ii], suffixes[jj]);

        // swap equal elements from right toward center
        uint64 numEqRight = index+length-1-pHigh; // number of elements equal to pivot
        uint64 numHigh = pHigh-pCurHigh; // number of elements greater than pivot
        if(numEqRight > numHigh)
            for(uint64 ii=pCurLow, jj=index+length-numHigh, numSwap=numHigh; numSwap>0; --numSwap, ++ii, ++jj) swap(suffixes[ii], suffixes[jj]);
        else
            for(uint64 ii=pCurLow, jj=index+length-numEqRight, numSwap=numEqRight; numSwap>0; --numSwap, ++ii, ++jj) swap(suffixes[ii], suffixes[jj]);

        uint64 numEqual = numEqLeft + numEqRight;

        if(numHigh > numEqual && numHigh > numLow) { // larger is largest partition
            if(numLow>1) multiQsort(T, dcs, suffixes, charIndex, index, numLow); // sort smaller
            if(numEqual>1) {
                if(charIndex<(dcs.v-1)/21) multiQsort(T, dcs, suffixes, charIndex+1, index+numLow, numEqual); // sort equal
                else qsortDCS(T, dcs, suffixes, index+numLow, numEqual);
            }
            if(numHigh>1) return multiQsort(T, dcs, suffixes, charIndex, index+length-numHigh, numHigh); // sort larger
        } else if(numLow > numEqual && numLow > numHigh) { // smaller is largest partition
            if(numHigh>1) multiQsort(T, dcs, suffixes, charIndex, index+length-numHigh, numHigh); // sort larger
            if(numEqual>1) {
                if(charIndex<(dcs.v-1)/21) multiQsort(T, dcs, suffixes, charIndex+1, index+numLow, numEqual); // sort equal
                else qsortDCS(T,dcs, suffixes, index+numLow, numEqual);
            }
            if(numLow>1) return multiQsort(T, dcs, suffixes, charIndex, index, numLow); // sort smaller
        } else { // equal is largest partition
            if(numHigh>1) multiQsort(T, dcs, suffixes, charIndex, index+length-numHigh, numHigh); // sort larger
            if(numLow>1) multiQsort(T, dcs, suffixes, charIndex, index, numLow); // sort smaller
            if(numEqual>1) {
                if(charIndex<(dcs.v-1)/21) return multiQsort(T, dcs, suffixes, charIndex+1, index+numLow, numEqual); // sort equal
                else return qsortDCS(T, dcs, suffixes, index+numLow, numEqual);
            }
        }
        return;
    }
    assert(false);
}

void BWT::qsortDCS(ntSeq const &T, DCS const &dcs, vector<uint64> &suffixes, uint64 index, uint64 length) const {
    assert(length>0); // qsortDCS should not be called with length 0 
    //if(length==0) return;
    // insertion sort serves as the base case
    if(length < 7) { // selection sort if partition has fewer than 7 elements 
        for(uint64 ii=index+1; ii<index+length; ++ii)
            for(uint64 jj=ii; jj>index && dcs.lt(suffixes[jj],suffixes[jj-1]); --jj)
                swap(suffixes[jj], suffixes[jj-1]);
        return;
    } else { // quicksort if partition has 7 or greater elements
        // choose pivot
        uint64 pivotIndex;
        if(length <= 50) { // if partition has 50 or fewer elements, select the median of first, middle and last elements
            uint64 pa=index, pb=index+(length>>1), pc=index+length-1;
            uint64 a = suffixes[pa], b = suffixes[pb], c = suffixes[pc];
            DCS::Compare comp(dcs);
            pivotIndex = comp(b,a) ? (comp(c,b) ? pb : (comp(c,a) ? pc : pa)) : (comp(c,a) ? pa : (comp(c,b) ? pc : pb));
        } else { // if partition has more than 50 elements, select median of 5 evenly spaced elements
            uint64 incr = length>>2;
            uint64 pivotIndexes[5] = {index, index+incr, index+incr+incr, index+incr+incr+incr, index+length-1};
            uint64 tmpIndex = findMedianIndex5(suffixes[index], suffixes[index+incr], suffixes[index+incr+incr], suffixes[index+incr+incr+incr], suffixes[index+length-1], DCS::Compare(dcs));
            pivotIndex = pivotIndexes[tmpIndex];
        }
        uint64 pivotSuffix = suffixes[pivotIndex];
        // swap pivot to the end
        swap(suffixes[pivotIndex], suffixes[index+length-1]);
        uint64 pLow=index, pHigh=index+length-2;
        while(true) {
            while(pLow <= pHigh && dcs.lte(suffixes[pLow],pivotSuffix)) ++pLow;
            while(pHigh > pLow && dcs.gt(suffixes[pHigh],pivotSuffix)) --pHigh;
            if(pLow>=pHigh) break;
            swap(suffixes[pLow++],suffixes[pHigh--]);
        }
        pHigh = pLow-1;
        // pLow points one beyond the smaller elements (the first larger element)
        // pHigh points one to the left of the larger elements (the last smaller element)
        swap(suffixes[pLow],suffixes[index+length-1]);
        // pLow now points at the equal element
        // pHigh points one to the left of the equal element

        uint64 numHigh = index+length-pLow;
        uint64 numLow = pLow - index;
        if(numHigh > numLow) { // larger is larger partition
            if(numLow>1) qsortDCS(T, dcs, suffixes, index, numLow); // sort smaller
            if(numHigh>1) return qsortDCS(T, dcs, suffixes, index+length-numHigh, numHigh); // sort larger
        } else { // smaller is larger partition
            if(numHigh>1) qsortDCS(T, dcs, suffixes, index+length-numHigh, numHigh); // sort larger
            if(numLow>1) return qsortDCS(T, dcs, suffixes, index, numLow); // sort smaller
        }
        return;
    }
    assert(false);
}

int BWT::verifyQsort(ntSeq const &T, vector<uint64> const &suffixes) const {
    cout << "Verifying qsort..." << endl;
    for(uint64 ii=1; ii<suffixes.size(); ++ii) {
        if(T.getSuffix(suffixes[ii-1]) > T.getSuffix(suffixes[ii])) {
            return 1;
        }
    }
    return 0;
}

inline uint64 BWT::getOcc(uint64 bwtIndex) {
    if(bwtIndex==term) {
        return 0;
    } else {
        uchar base = bwt.get(bwtIndex);
        uint64 occIndex = bwtIndex/INDEX_PERIOD;
        uint64 tmpOcc = occ[base][occIndex];
        uint64 startIndex = occIndex * INDEX_PERIOD;
        for(uint64 ii=startIndex; ii<startIndex+bwtIndex%INDEX_PERIOD; ++ii) {
            if(bwt.get(ii)==base && ii!=term) ++tmpOcc;
        }
        return tmpOcc;
    }
}

inline uint64 BWT::stepLeft(uint64 bwtIndex) {
    if(bwtIndex == term)
        return 0;
    else
        return C[bwt.get(bwtIndex)]+getOcc(bwtIndex);
}

vector<uchar> BWT::unpermute() {
    vector<uchar> tmpSeq(n);
    uint64 nextIndex = stepLeft(term);
    for(uint64 ii=0; ii<n; ++ii) {
        tmpSeq[ii] = bwt.get(nextIndex);
        nextIndex = stepLeft(nextIndex);
    }
    reverse(tmpSeq.begin(), tmpSeq.end());
    return tmpSeq;
}

int main(int argc, char **argv) {
    //ntSeq T("ACGTACGTTTGGGGGTTATATAAATATAAATTT");
    //ntSeq T("ACGTTACCG");
    ntSeq T("ACGTACGTTTGGGGGTTATATAAACCCACACCCCTTATATTTAAAATATATATATTTTTAATATAATATAACCCCCCATTCCCCGGGAGATTGAGTCCCCACAAT");
    //ntSeq T("ACGTACGTTTGGGGGTTATATAAACCCACACCCCTTATATTTAAAATA");
    //ntSeq T("AAAAAAAT");
    //ntSeq T("AAAAA");
    cout << "T.size(): " << T.size() << endl;
    cout << "main T: " << T.getSuffix(0) << endl;

    uint64 numBlocks[] = {1,2,4,9,10};
    for(uint64 ii=1; ii<10; ++ii) {
        for(uint64 jj=0; jj<4; ++jj) {
            uint64 logv = ii;
            BWT(T, logv, numBlocks[jj]);
        }
    }
}
